!> 预估校正法(半隐式)比龙格库塔法更适合刚性常微分方程求解
program example_stiff

    use seakeeping_module, only: runge_kutta, predictor_corrector, wp => sk_real_kind
    implicit none
    real(wp) :: t, dt, x(2)

    t = 0
    x = [2, 0]
    dt = 0.005_wp

    print *, "****** Predictor-Corrector ****"
    print "(2(A10,1x))", "t", "y"
    ! 预测-校正法
    do while (t < 0.05_wp)
        call predictor_corrector(fcn, t, x, dt, 1, 2)
        write (*, "(ES10.3,',',ES10.3)") t, x(1)
    end do

    print *, ""
    print *, "****** Runge-Kutta ****"
    print "(2(A10,1X))", "t", "y"
    t = 0
    x = [2, 0]
    dt = 0.005_wp
    ! 龙格-库塔法
    do while (t < 0.05_wp)
        call runge_kutta(fcn, t, x, dt, 1, 2)
        write (*, "(ES10.3,',',ES10.3)") t, x(1)
    end do

contains

    subroutine fcn(t, y, yp)
        real(kind=wp), intent(in) :: t      !! 自变量
        real(kind=wp), intent(in) :: y(:)   !! 应变量
        real(kind=wp), intent(out) :: yp(:) !! 导数

        yp(1) = y(2)
        yp(2) = 1000.0_wp*(1.0_wp - y(1)*y(1))*y(2) - y(1)

    end subroutine fcn

end program example_stiff
